# Figure5_Assessment of drought damage and recovery.R
# Coded by Prof. Alastair J. Potts
# 31 Jan 2025
##########################################
# Use groundhog package for reproducibility
use.groundhog <- TRUE

if (use.groundhog) {
  library(groundhog)  # Load groundhog package
  set.groundhog.folder("C:\\Groundhog_R\\R4.4.1_2024-05-14")  # Set groundhog folder location
  GroundhogDay <- '2024-05-14'  # Define GroundhogDay version
  groundhog.library("tidyverse",GroundhogDay)  # Load required packages using groundhog
  groundhog.library("lubridate",GroundhogDay)
  groundhog.library("drought",GroundhogDay)
  groundhog.library("patchwork",GroundhogDay)
} else {
  library("tidyverse")  # Load required libraries normally if not using groundhog
  library("lubridate")
  library("drought")
  library("patchwork")
}


# Set working directory
setwd("D:\\Dropbox\\100_PROJECTS\\2024_ThicketHydraulics_SkeltonButtnerPotts")

read.csv("1_data/08_InFieldCanopyAssessments.csv")->
  dat

dat%>%
  mutate(FoliarCover=100-Defoliation)%>%
  filter(Species!="Boscia_oleoides")-> #Boscia oleiodes was in a leaf flush and could not be assessed with the metrics 
  DAT

#### MAIN FIGURE ######################################################
DAT%>%
  group_by(Species)%>%
  summarise(p_value=t.test(unlist(Crown_extent)~unlist(Year))$p.value)->
  ce_pval
DAT%>%
  group_by(Species)%>%
  summarise(p_value=t.test(unlist(Defoliation)~unlist(Year))$p.value)->
  fc_pval
DAT%>%
  group_by(Species)%>%
  summarise(p_value=t.test(unlist(Dead_leaves)~unlist(Year))$p.value)->
  dl_pval

quickplot <- function(obj) {
  DAT%>%
    filter(Crown_extent!=0)%>% # think of how to deal with the dead plants...
    rename(VAR=obj)%>%
    mutate(Species=gsub("_","\n",Species))%>%
    #mutate(Species=factor(Species,levels=c("Poyul")))
    
    group_by(Species,Year)%>%
    summarise(uCE=mean(VAR),
              seCE=sqrt(var(VAR)/(n())),
              sdCE=sd(VAR))%>%
    mutate(semin=uCE-seCE,
           semax=uCE+seCE,
           sdmin=uCE-sdCE,
           sdmax=uCE+sdCE)%>%
    ggplot(aes(x=paste(Species,"\n",Year),y=uCE))+
    geom_errorbar(aes(ymin=sdmin,ymax=sdmax),size=1,width=0,colour="darkgrey")+
    geom_errorbar(aes(ymin=semin,ymax=semax),size=2,width=0)+
    geom_point(size=4)+
    geom_line(aes(group=Species))+
    theme_bw()+
    xlab("")+
    geom_vline(xintercept = seq(2.5,by=2,length.out=4))
}

quickplot("Crown_extent")+
  # annotate("text",x=seq(1.5,9.5,by=2),y=rep(100,5),
  #          label=case_when(
  #            ce_pval$p_value<=0.001 ~ "**",
  #            ce_pval$p_value<=0.01 ~ "**",
  #            ce_pval$p_value<=0.05 ~ "*",
  #            ce_pval$p_value>0.05 ~ "ns"
  #          ))->A
  annotate("text",x=seq(1.5,9.5,by=2),y=rep(95,5),
           label=round(ce_pval$p_value,3))->
  A

quickplot("Defoliation")+
  # annotate("text",x=seq(1.5,9.5,by=2),y=rep(100,5),
  #          label=case_when(
  #            fc_pval$p_value<=0.001 ~ "**",
  #            fc_pval$p_value<=0.01 ~ "**",
  #            fc_pval$p_value<=0.05 ~ "*",
  #            fc_pval$p_value>0.05 ~ "ns"
  #          ))->B
  annotate("text",x=seq(1.5,9.5,by=2),y=rep(95,5),
           label=round(fc_pval$p_value,3))->
  B

DAT%>%
  group_by(Species)%>%
  summarise(
    p_value=t.test(unlist(Leaf_discolouration)~unlist(Year))$p.value)->
  ld_pval
quickplot("Leaf_discolouration")+
  annotate("text",x=seq(1.5,9.5,by=2),y=rep(80,5),
           label=round(ld_pval$p_value,3))->
  C

quickplot("Dead_leaves")+
  annotate("text",x=seq(1.5,9.5,by=2),y=rep(50,5),
           label=round(dl_pval$p_value,3))->
  D


# SPI Figure
qCols <- rev(c("#3B9AB2", "#78B7C5", "#74A089", "#E1AF00", "#F21A00"))
read.csv("1_data/03_CHIRPS.csv",skip=3)%>%
  as_tibble()%>%select(-X)%>%
  mutate(ym=ceiling_date(ym(paste(Year,Month)),"month")-1)%>%
  ggplot(aes(ym,spi_val))+
  geom_vline(xintercept=ymd("2020-SEP-30"),size=2)+
  geom_vline(xintercept=ymd("2022-AUG-25"),size=2)+
  geom_segment(aes(x=ym,y=0,xend=ym,yend=spi_val,color=spi_val),size=2.5)+
  scale_color_gradient2(midpoint=0,high="blue",mid=("green"),low=("red"))+
  geom_line(size=0.5)+
  theme_bw()+
  scale_y_continuous(breaks=seq(-3,3),limits=c(-2,2))+
  scale_x_date(date_breaks = "year",date_labels = "%Y",
               limits = c(ymd("2016-01-01"), ymd("2022-12-31")),
               expand=c(0,0))+
  ylab("12-month SPI")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.position="none")+
  xlab("")+
  geom_hline(yintercept = 0)->
  Fig0


Fig0/
  (A+ylab("Crown extent")+scale_x_discrete(labels = NULL, breaks = NULL))/
  (B+ylab("Foliage (%)")+scale_x_discrete(labels = NULL, breaks = NULL))/
  (D+ylab("Dead leaves (%)")+scale_x_discrete(labels = NULL, breaks = NULL))/
  (C+ylab("Leaf discolouration (%)"))+
  plot_layout(heights=c(1.2,2.2,2.2,2.2,2.2))

ggsave("3_results/Figure5.jpg",units="cm",height=20,width=20)
